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Abstract 

W avelet bases are used to generate spaces of approximation for the resolution 
of bidiiueosional elliptic ini parabolic problems. I nder some specific hypothe- 
ses relating the properties of the wavelets to the order of the involved operators, 
it is shown that an approximate solution can be built. Tins approximation is 
then stable and converges towards the exact solution. It is designed such that 
fast algorithms involving biorthogonal multi resolution analyses can be used to 
resolve the corresponding numerical problems 

Detailed algorithms are provided as well as the results of numerical test? cn 
partial differential equations defined on the bidunensional torus 
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Introduction 


Variational approximation methods for partial differential equations are based 
on weak formulations and on suitable spaces of approximation. 

Wavelets are known to he unconditional bases for a large variety of spaces 
and therefore are good candidates for the generation of approximation spaces 
for partial differential equation problems. The goal of this paper is to show that 
moreover, wavelet bases may lead to fast and adaptive numerkal resolution of 
tbe corresponding approximations. 

Iq this paper, as in previous papers (J. Liandrat and P. Tchamitchian [10], 
(11'). the wavelets are used to expand the approximated solution of a partial 
differential equation as well as to approximate tbe kernel of the differential 
operator. They are not used only to perform “the linear algebra' (O. Bey 1km 
(2)) related to more classical methods of resolution. 

Starting with an expansion of the form / = £*</♦;,>**. the solution 
of the equation Lv - f where L ts a constant coefficient elliptic differential 
operator is u(x)s //(y)A'(r.y)<fy where A*(r.y) = 

I nder suitable conditions that will be made precise later, the functions 
!**♦» as well as L'tx are pseudo wavelets, very clone to wavelets (V Meyer 
1 13] i. This turns out to provide a stable approximation of «. However, the effi- 
ciency of the corresponding numerical approximation of v relics, at least in this 
work, on the hierarchic structure of niulUresohition analysis since it provides 
fast tree algorithms. We will shew that, if the operator satisfies suitable condi- 
tions that will be made explicit later, then the above mentioned pseudo wavelets 
are directly related to biorthogonal multiresolution and wavelets I'nder these 
conditions, competitive numerical algorithms involving 0(.Y t or Ot.Y log .V) op- 
erations can then be derived. 

This paper provides the analysis of the problem and exhibits the '■orrespo.id- 
ing numerical scheme* It is then organised as follows. 

The first part is devoted to the general concept of hiorthogornl multitesoln- 
tion analysis on £ J ( J?’’). In the second part focus on the problem of the sta- 
bility of the multiresoluuon framework under the action of constant coefficient 
elliptic operators. The cases of homogeneous and non homogeneous operator* 
are tteated separately The third part deals with the numerical algorithms while 
the last section is devoted to numerical tests related to the resolution of elliptic 
and parabolic equations in bidimensional spaces 


I Generalities: Biorthogonal Multi-Resolution 
Analysis in L^JR") 

The concept of multiresolution is at the basis of our construction and we there- 
fore start with a short description of it : 
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Definition 1.1 (Y. Meyer. (12]> 

A T*rxf%ltr mntUre solution analysts of l 3 (ST) is a sentence of increasing 
cicftd smispaets I}, j € Z. Ij C Vj-i satisfying Ike follounng conditions 

»' f£* '!> = {*>)• U!£ V » d<*«* « L 7 [tT): 

»*/ /(*) € f(Zr) € V; +l 

««.* f(x) e l‘o «* /(x - k) 6 r : . V*- € Z’-: 

*Arrr <rwf$ a function ♦ ia t' 0 . svrA <*■/ lAc s«f of function* {♦(x-Ah k € 
Z”). is • Jtiesr busif * /or to: 

r.i tke function ♦ tj regular and localised : ♦ is C \ ♦ is almost 
titryvhtn iiffertntiaUe. and for almost trtry x £ ft* . for fiery integer 
o < r and fcr all integer p. i { tints C t such tkal 

< C,(l + IxiH (!) 

A consequence of u). in) and trj is that each 1} is generated by the family 
ot‘ functions {♦.*(x) = -It), A € Z*| 

For simplicity reasons, wr will only consider the case a = 2. but all the 
results presented to this article can be generalised in any dimension We trill 
always for vectors a contracted notation if e is a bidiinensional vector then 
e = 

1.1 Orthogonal multiresolution analysis 

To build an orthonotmal muitiresolution analysis, the Ricsr basis {♦< -*'. A € 
Z'} is first otthonormalixed in such a way that the resulting etthononnal basis 
is still of the form {♦;. - k).k€ Z 3 }. 

The wavelets are introduced via the orthogonal complement of 1> in \\ : ||'-. 
Mcee precisely, if E is the set of all vertices in the unit cut* [C 11 v and if 
r* = r \ {0} . we have the following theorem: 

Theorem 1. 1 Tktrr an X functions 'cf is H’ 0 such that the collection 
A € X 3 .! 6 E ' } is an orthonortnal iasis of lip. Monootr. tack ♦' 
satisfies the same property (If of regularity and localisation as ♦ and moreoter. 
satisfies the following cancelation property 

3 m € £V . such that V A = (Aj A 3 > 6 X 3 . 0 < A, < m, 

( * <?' 

/ r'Vlrldi = 0 

Jm ! 

'A roiv.sion of vectors {f j \ € A). in • IHH-eri »|>»ce H m a R»e«» h*j,« if T-r-cr 
r € ff can be nritien in a urjuiie may w i * IT " l d «h*rv (£ *«ai 3 }*•'* i« Unite »r*d 
defines a r^rm •‘jniraJe.u u, ;|x 
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The scaling invariance property n) implies that, for til >, the family k £ 
J?\ f 6 f} is an orthooormal basis of IV;. We will abo use the following 
contracted notation: , A € \i\ where A. = {A a 2 *' (t + 5) . * € Z*,f € 

£ * } Indeed, there is a straightforward Injection between A; and the set of pairs 
{(£.*).* 6 Z 3 ,s € £*} We will also use the following sets: A « Uj«»A; and 
AI-JT..A,. 

From the inclusion l’c C l'i the following set/tag re/afto* can be derived: 
+(r) s £ h M* ~ l) • 

while from U”? C li one obtains the following detail rc/tfroa* 

#'(T) U Y ^♦(2x - /). Vt € r* 

»«*» 

!r very useful to transform these relations using the Tourier transform 
*hicli is gi\Yn by the equality 



Indeed, the scaling relation* then become 

*l« = Vc«/2)^/2) (3) 

and 

*MO= M,^/2v4«/2) (4) 

ahere Afo(£) = h’C* 1 f ' and A/ f (^ = dL;**' C°° 2r peri- 

odic functions. 

This leads to 


1*1 <*• * 


I.V 


The following conditions arc satisfied 


Y '* + "1 = • V (r.r'l 6 £ 3 (6' 

•ft 

and 

A/,(Tei = A„ V(f.ft6£ J <:» 

and nre called, following electrical engineering terms the fuadnifm miTtwrfihtr 
conditions \\e will also call the functions M, quadrature rtnrccr filters 
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Conversely, it is shown in A. Cohen et al [4j that four 2r periodic functions 
A/,( 0> f € E satisfying the quadrature mirror filter conditions (6) generate 
through (o) an orthogonal multiresolution analysis if some specific conditions 
are satisfied. 

Remark: 

• In this paper, we will often use specific multi-resolution analyses of l 3 (JR 1 ) 
based on a tensorial product of multiresolution analyses of More 

precisely, such analyses are defined as follows: if (V ; ) is the sequence of 
spaces of a ID multiresolution analysis and if >V ; . y, v. mo and n» t are re- 
spectively the related wavelet spaces, the scaling function, the associated 
wavelet and the quadrature mirror filters, then, the sequence of spaces 
(1} k. defined as 1} = Vy ® is a multiresolution analysis in J? 5 . More- 
over. ♦-» = i?fr ; Vix?) is the corresponding scaling function; (IT, ). with 
li j — 1 1;* ate the wavelet spaces; the three generating wavelets are 

*o.(x) = *1x,)v<x 3 ). *; ,(*)=s v(* : )v(r 3 ) and : (x) = i*(xitax 3 i; 
— rn t i l ift )n r, « J ff 3 ) with f € £ are the quadrature mirror fillers 

1.2 Biorthogonal Approach 

A relaxation of some properties of orthogonal multiresolution analysis can be 
performed using the one biorthogonal approach. This approach provides some 
flexibility since it allows to distribute the relevant properties of the mult ires 
oiution (number of ter© moments, compact support or regularity) to the two 
involved multiresolution analyses. Moreover, it will turn out to be that the 
biorthogonal framework is “stable'’ under the action of a large class of opera- 
tors while the orthogonal framework is “fragile" . 

Definition L2 tfc nil ltcr1keyp»*J nnhm solution analysts of L ! {Mt 3 ). two 
mahtrtsolultop analysis (( ;);$* (f j lid sacA tAat thtrt enststvo famihts 
of corrtsfoudin^ scatiwf functions r and ? inch that r. >.) s ( it . for 

oil j j'e 2 T and. k and *’ € M 3 . ' 

In this ciue^we define the wavelets spaces A, and X } a* (\+: = l } + A'.. 
= fy ^ Xj with L XX, 4-J.A'j aiKl we introduce the functions (*\ = 
Vfi t {2 } . - k) and fix — 2 J fi't’P — A), ; € f* that generate respectively .V, and 
X. and such that (fix. ?»•) s ( t i-( }} <fkk- 

Moreover. following the construction of orthogonal multircsoluticti analysis 
we define 2x4 filter* (i e. C°° 2» periodic functions). P, and P,. t € E. 
associated with the two biorthogonal multiresolution analyses. These filters 
satisfy the biorthogonal quadrature mirror filter rHations equivalent to (fii that 
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are 


V{.{‘ and ’0, 2»j* , 

£ •+ T ')*Vtf •+ «) = . 

tit 


it) 


^»i *«^r> »> • 


As in the orthooormal case, the generaliratioru of relations (3) and (4) relate 
the scaling functions and the wavelets to the filters as 


= fttf/2)?tf/2) * n*™ W 




?tf) 


= PoW'Msu) = 

j=> 


■sf 

* tf) = 


ntf/21?tf/2> (101 


As in section 11. the question to know is under which conditions the biort hog- 
onal filters P, and P t satisfying the quadtature mirror filter conditions (>1 de- 
fine two biorthogooal multi resolution analyses’ Again, a specific condition is 
required and has been formulated in A. Cohen et nl. A\. We will use a weaker 
version of this formulation adapted to the case of functions with fast decay. It 
can be expressed in term of the following theorem for which a complete proof 
can be found in Pj. 1‘onenti [161 and \V. Dahmen and A Michelli (5) 


Theorem 1.2 lei a > 0 andle and P,{() ( € E.k cijhfC’ 2* periodic 
functions satisfy mg tkc kiortkogonal yuadratorr mirror filter* condition * ; S). 
Defining r, r, 0, and $ using formula (9) and (10 j. i / 


- tkrrt enst C and jj > 0 sac* Ik at for all ( £ JR 3 

l ? tf)i < c(i + *ir w . < ni + £jr'-* . (in 


- tr *.<#**€ Z 3 . 


• ini i/ 


thtn 


\r{x - k).r[r - k’)) = f ki .. 


/«> K: J ' iCF * 


< X 


1121 

(131 


• tke segue uos of euispace* ((•) ) ;ig . and (f? ) j6 * generated respectively 

b {”•* • * € X" 3 ) and 6 X J ) art two kiortkogonal mthirt solution 

in<tly$<$: 


* Ik e n »’ e let families { 4 x (r) = V*'{Vt - *). X € A > and = 

Ve , (Vr~k':. A € A} orr in o kiortkogonal Ries: bases o/L 3 {JR : ). 
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II Constant Coefficient Elliptic Differential Op* 
erators and Biorthogonal Multiresolution Anal* 

ysis 

II. 1 General results 

The starting point of this section is the following remark. Given (♦*) a family 
of orthonormal wavelets and knowing / = £{/. the solution of the 

equation 

/.us/. (141 

where I is an elliptic operator of order t is. at least formally. 

!♦*!■ i 151 

When I and L~' are bounded on i 5 {iR J ). the families { t ~ 1 ♦ x } *nd 
ao* two biorthogonal Riesr bases 

The question we address now is related to what happens in the specific case 
of wavelets when the operator L is unfunded (aa in the case of a differential 
operator). 

la the following paragraphs, we first see that, assuming some compatibil- 
ity conditions between * r and i. the two families of functions and 

} are wavelets or pseudo wavelets (Y. Meyer '13' > depending on whether 
the operator is homogeneous or not Then, we show that in some caaes. a 
biorthogonal framework embedding {/.“’•♦*} and {!*•»} can be built 

Tobe more precise, we take Hi) an t-regular niuiiiresolui’.on of con- 

structed using a leitaorial product of two ID tnulUresolutions and L a constant 
coefficient elliptic differential operator. Let us write L - £ l<a <, where 
D is the operator We define in a standard wav the symbol of I as 

<-!() = (16) 

It is a polynomial in £ of degree s and. if / 6 C'ff? 7 ), we have the well known 
formula 

Jf 5 .Z/«) = /UW <17i 

We no it formally 

P T rr and£ = L~'r 

and more generally 

#\ = 2 and^ = ¥ 9 L- % *\ ,19;. 

Note that when L is homogeneous - i) and (x) s 

¥P{V* — t) *hJe. iu general this is not uue. 
i hen we have 



Theorem 11.1 G»re» « ftmi/y efr reftler uoreUts u-tf* m + 1 zero moment*, 
if L ij t homoftntfius operator or if J. is <i inhomogeneous open tor mtk o 
stnctlf positive spmkot of order $ > 0 ond sock tkot r > s ond m > s. tken 
RctuUntr r € ir- and ff 1 € /T" 

/ ocehnticn for oil multi- indict r >' and ■> sick that jr' < p - * and 
hi< r + * - 1 and oil integers I € A*. 

I |^'> x (x)i < c, 2»< v, y (i + y|r-A|)“ ; 

\ |fl'#l(x)| < C, (l-r>- A|)’ 3+ '”"" ' 


Ceneeletton. Let s* s *{»). 

x*0*(xWr = 0. 0 < 1*1 < n> 4 a 

t'frtrUr = 0. 0 < (i| < n> — « . 


Let x* : 

iii 


there assn the komoftn eon* cast end a — 0 in the inhomogeneous one. 

A complete proof of this theorem can be found in V. Meyer [13] and Ph. 
Tchamitchian (17], 

Remark: 


• Following A Meyer [13] and according to Theorem 11.1. a factorization 
of the operators i* and L~ l can be performed as L’ = f o T” and 
= Co T* where T is a diagonal operator in the wavelet basis defined 
as: r v •- Vi* aad where C and C are bounded on i. 3 and defined by 


f’ : Va*-*$a . C : «— • (W) 

The operators C and C act just as a traiudormation between two bases 
1 he operator T is nothing else but the classical preconditioning operator 
for elliptic problems (S. Jtfttd [£jl that mimics a diagonal deritation in 
the wavelet basis 

In other words, thanks to this factorisation, the computation of the image 
of a function by an elliptic operator or its inverse can be transformed into 
a well conditioned problem using a diagonal operator in a suitable wavelet 
basis. This is essential since it provides the numerical stability cf tbe 
further developed algorithm. 


We can then rewrite 1 13) as 

U = Z 2 ~ :, 'f 



* We call a mwki-iiwtev snv couple of iweger v * 


( 21 ) 



An important issue, as far as numerical applications are concerned, is the 
computation of tbe sum (21). Indeed, even if this sum corresponds to a pseudo 
wavelet decomposition, fist algorithms for the computation of the sum are not 
available, la the framework of this paper, the fast algorithms are linked to 
the concept of muhireaolutioo analysis presented in section I We prove in 
the following sections that, under suitable conditions, the construction of a 
multiresolution analysis embedding the function t*i is possible. Moreover we 
provide explicit expressions of the quadrature mirror filters required for the fast 
implementation of (21). 

The starting point of our construction is due to P.G. Lemarie [9] who con- 
structed in the one dimensional case two biorthogonal multiresolution analyse* 
from an original orthogonal one and from the derivata operator. We generalise 
this approach to any dimensions and for any homogeneous elliptic operator. 

In contrast to the classical constructions of multnesoluiion analysis, this is 
an inverse problem. Indeed, knowing the two dual wavelet bases we car. define 
two sequences of subspaee* ( A*. ) and (A. ) 

A; = span(dx. A € A,./ < ;} 

A. = $p*n{ . A € A ' 

The open question is the following: how ran we construct two sequences of sub- 
spaces (( \ ) and (f, ) for which ( Xt 1 and (A, 1 play the role of two bier t hogonal 
wavelet spaces'* 

In other words, the problem isjhe construction of the generalised scaliug 
functions r and t related to 0 uid 6. 

In the case of non homogeneous elliptic operators tbe approach used in the 
homogeneous case can not be transposed and we will not be abie to define a 
multiresolution framework embedding tbe space sequences (£", ) and (C\ V How- 
ever. we will show that the essential property of embedding spaces as well as 
the existence of scaling (31 and detail (4) relations can be saved. This will allow 
us to derive fast and stable algorithms to sum up (21) even in the case of non 
homogeneous operators. 

II.2 The case of homogeneous elliptic operator 

In that case the natural candidates for r: !"•♦*. are not defined in L' for the 
basic reason that • does not belong to the range of L. or. in other words, suffers 
from a lark of zero moments Using a preconditioning operator, we will adapt 
the function 0 to the operator L~ l (i.e. we will transform ♦ so that it enters 
the range of i) while preserving the two scale relations (3) and (4) We will 
then check that the resulting multiresoiutioo Analysis is fitted to the functions 

and previously defined 

More precisely, we have the following theorem 



Theorem 0.2 Lti ( Vj ) be tie fomiip of tmiedded spaces of «* r-rtfnler nut- 
tmsoUtion analysis ofl\ I?*), lti L it t komopencons operator of orders end 
of symbol c an d lti S it a 2* periodic function, not tontskinp on ]0.2»p end 
tfonoUnt to o tn uro. Ifm > a- 1. r > 1. ond if Ike etpki functions />,({) 

nod P,($). f € T defined os 


Pc (0 ■ 2* 


5(g) 

5(20 




= 57 


1 5(2* ) 


5(0 


•'MM . 


(23) 


P,(0 = 2’ 5(0 U,(0. A(M = (24) 

2* 5«) 

aft C*. a > 0 then they satisfy tie ynadrainre mirror fitter conditions (8) end 
define tioo htorikopennlmnHirtsolthon analysts. Tkt corrtspondinf novelets 
ore tkt functions 6 ond f end tkt svaltnp functions r end r err dent ed follomny 

(*/ and (Kfl 


Proof: 

O By construction all the filters a re C* with some a > 0 and they satisfy 
the hiorthogona) quadrature mirror filter conditions (d) 

The only point to prow is the convergence in i J (i? 2 } of the infinite prod- 
ucts (9) sad (10) defining the two scaling functions We will use the following 
lemma: 


Lamina II.1 Iff p(x). x € W 3 . it a homoftntons polynomial of degree s. 
ond let S end C, it 2r periodic functions : tien tkt follotcmp propositions on 
rf»r>-«/« nt 


it 




5(x) 

P(r) 


(25) 


ni 

5(x) = 2*Ctx/2)5(x/2) . (26) 

5< x>,~ : p(x) (27) 

Proof: 

O The equation (26) is obtained front (25) wntten for x and x/2. while (27) 
is derived from (25) when x — 0 since necessarily C(0) = 1 
Conversely. (25) is obtained from (26) Indeed, since 


5lx) 

FT) = Cu/2, FT72) 


5(2- v x) 

p<2- v x) 


aV 


nw- 


X? . 


L 


thanks to (27) are obtain (25) when X — oc ■ 


This lemma allows us to calculate the infinite product (9) and (10). and to 

ho = n ***"'« = 

= n ^2**o = n f(2*-’on^(r , (iK(2-''( J i . 


aud finally. 


?«* = ♦*> ■ no . || ♦«. 


Ihe function ♦ being r-regular the conditions (11) are trivially satisfied. 
Defining the function 0 et 0 by (9) and (10) we get 

= ^((/2)?U/2) = 

5« .. if.-i (29 j 

* (0 = i*e«/2)?«/2) = . 

The wavelet admissible condition (13) is immediately satisfied thanks to Theo 
rem II 1 Finally the assumption (121 is satisfied by construction, that completes 
tht* proof m 


Remarks: 

• Note that here, thanks to homogeneity, the subscript * recovers its classical 
“wavelet” meaning since in thai ease. <**!») = 2 i/J *(2-' i x-b) and fi x (x> = 
2> ! -9{{2-'x-kl 

• The relation (25 'i is a general nation in any dimension, of the classical 
formula 

ft = SSlI 
•*1 * 

used by P.G. Lemarie for the first order differential operator in one di- 
mension. 

• The function 5 can be interpreted as the symbol of a difference operator 
that we will call Dl If 5 is a trigonometric polynomial, then D t is a 
finite-difference operator and 5 is C The fact that 5(x) - e{r) is 

just the translation that Dt is consistent with L Moreover * S remove* 

exactly the singularity of t for ( s 0. Conversely for r. the singularity 
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given by 5 at points (2tn.2:rn). n € Z. will be removed exactly by o in 
mo and by the teros of ♦ at points (2rn. 2ml, n € X ' . Notice that this 
last point won't be true for inhomogeneous operators. 

From a certain point of view, Di can he seen as a preconditioner for L 
since r and ? are defined by 

r= i-'tft^andf *(0£r ; r*. (30) 

• In one dimension there is a canonical choice for 5 and therefore for Di 
such that, if the function ♦ and • have a compact support, then r. ?. 
(?, and 0 are also compactly supported. Indeed, in that case we have 
necessarily o({) = o£*,a € € Therefore, the canonical choice for 5 
is 5(4) = «(-«)’ (1 and Di is then a non-centered finite- 

difference approximation of £ of order I. Indeed, it is well known (see 
I. Daubechies (6]) that the quadrature minor filters related to orthogonal 
compactly supported functio ns 4 and ♦ are: m<)(0 = (1 4 e'') m £(£) 

and mi(4) = e ,f (l - e ~ s( ) £(£ -*• -) where £ is a finite trigonometric 
polynomial Then we get 

*><0 s 2* (HV< f '*£(£). 

(l-«* )‘C[*+T). 

> 0 (O = 2- (I + ««)" (1 + *-<)' £(4). 

(1 - + T) 

which proves that Pv Pi. P® and P\ are also finite trigonometric polyno- 
mials Then, using the following lemma borrowed from G. Deslauriers and 
S. Dubuc [7] we deduce that the functions ?. r, 9 and 9 have compart 
support. 

Lemma n.2 

V no = E, v *.v, «*'* E-i*. v « 1. thfn n~ , r(?-^) is 

entire fnntUo* of eipontnU*! fjj*. !n ptrtuuUr it is tht Fenner Irmas- 
form of t itfilnbtUo* ri ih swppori is [A*j. .Vj) 

Clearly. this canonical form is no longer available if the space dimension 
is larger than 1 since the multidimensional quadrature mirror filters can 
not be factorized art above 

II.3 The case of inhomogeneous elliptic operator 

Here the non homogeneous property of the operator is obviously not adapted 
to the scale invariance property of the multiresohition analysis. We will see 
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however, that introducing at each level a new scaling function, an embedded 
family of spaces can be constructed which preserves the mathematical properties 
relevant for numerical applications. 

The natural candidates for r* are a the functions £"'♦» They arc well 
debited iu but suffer now from a lack of locaJiratton when j incream 

Indeed, we have 

Jim^ i iv ; i] (i*) - G(/ - 2~*i) i = 0 

*?” — t, 

where C is the Greens function of the operator L defined as G{{ ) = 1 .• i fc*r 
example, when l - 1-A G(0 = l/U+^laod Girls c* > C dccream fast 
but mathematical and numerical difficulties come from the fact that the family 
of functions { G\t - k'2~ ! ).k 6 Z 5 } is not a set of function* rescaled with j (in 
other words, this family is not obtained by rescaling nod translation of a single 
initial function! This implies that the control of the localisation by the index 
j «s lost It follows that the family {L~ * x . t e Z 3 } » not a g w J basis to 
reconstruct our solution 

Let us show now that, however, a process very close to the one used in 
the homogeneous case will provide an efficient algorithm, for the summation of 
formula (21) 

We mimic the construction performed in the homogeneous ease Let . I' . 
be an r-regular xnultirt solution analysis let L be an elliptic operator of orders 
with constant coefficients and <r its symbol (we now suppose <r({ ) > > (i ye) 

let us also define the homogeneous polynomial of order s. v as the principal 

part of V. and let 5(f) be a 2r periodic function with >’;{ ! »-{'•. where 

J — c 

n will be fixed later 

Then. Vj £ Z. we define a difference operator D, by its symbol >(f . I 

Following the previous section, we define Vj. k € M 

?.* = y'l-'D,*., til: 

and 


i32i 


= y/'M*,*. 

By construction, and thanks to the fart that / is a constant coefficient 
operator we have r »{r) s f } \r - and f.»(ct = 0, ,r - £) where 


and 


<r,( 

- , . iic 2 ' ■ 

p.io = y~' — — 


(.IT* 
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Remark 


• The functions tj mimic the function ? defined in (28). liifortunatelv. it is 
not possible to define the equivalent of r (28) since D' 1 LQ;k € i J . Note 
however thnt. by chance, (lot involve* directly the Sx functions 

Then, with Pt- and P t defined in (23) and 24 we get the following scaling 
and detail relations: 


W = 2 )?._.(?) 

£($) = 2A«/y*‘)?, + ,«) V f €£' 


(35) 


Remark: 

• An important point is that the fillers P, are independent of the scale 
index j as it is originally the case for standard multiresolution analysts. 
Furthermore . since they are defined by (23) and 24. the filters P, are 
directly related to the homogeneous operator L cf symbol & if n * s 
This poirt is essential sic.ce it means that if D_ is consistent to L. then 
the tree algorithms related to the niultiresclulion spaces ( f y ) and used to 
sum up (21) are stable even if the functions e. ( r ) are not standard scaling 
functions. Indeed. 5(f)Mfl — <_.j 0 and therefore ?j(0) = 0 

In other words, even if the functions ?ys(x5 are used as scaling functions on 
the range of L. they haw rero moments as wavelets have 

Finally, we can prow the following theorem. 

Theorem 11.3 Fcr 0 < n <*••,#< r aai s < m Iht /aeefion* ?,> defined 
tf (S3) tad (31) satufy 


!£’?..<*) < O, («*?.> - {36) 


If n — s and i if £ £ £’ art ( ° then f table tret algorithm* are ai^ailiible 

re. 

30 < C < C < +x itch that if f - 3\0\. then 


C <||/ji : < C £|d* 2 . 


(37) 


Proof. 

• Since 5 e C* the we can apply theorem 11 1 that proves the localisation 
inequality (36). 


• The relations (36) defines directly the tree algorithm? required to compute 
the coefficients {e ; *} such that 

and since, when » s s, the involved quadrature mirror filters could be 
related to the homogeneous operator L the stability of the algorithm is 
equivalent to the stability of the transform {c;*} — / for 

• This transform is stable if and only if the family {?«*. t e E 7 } is a Riesz 
basis of A’?. 

We have 

ii/r = 


Since {♦,*} 1 $ a Riesa basis, since 5(0 i* bounded, and since e is bounded 
from below. ihen ! : / H 1 ^ which is the second part of the 

inequality 1 37 ) 

To prove the first part of (3?) we use again the fact that the filters P and 
P, are related to L 

Indeed, if we define 4* replacing o by & in ( 29) and if we define / as 
f — ]T * d\0\ then the transform / — » {d* } is stable Moreover thanks to 

theorem 11.1. the operator / — / which can be also defined as *> — 
is also bounded (Y. Meyer [13]) Therefore the operator / — (J, } is 
bounded that is the first pan of (37) and that completes the proof. 



Ill Approximation and Numerical Resolution 
of Elliptic Problem on the Torus 

Thi* action is devoted to Uk approximation of elliptic problem* on a sequence 
of embedded Galerkin spam awooaled with a multiresoluticn analysis anJ to 
tht corresponding numerical algorithm* 


1 ! 



Classically, w * will consider the problem with periodic boundary conditions 
to avoid the difficulties of general boundary conditions. We will use a r regular 
multiresolution analysis of the torus 7j* jp = (ft/M)* constructed using a clas- 
sical periodiiatioc technique (Y. Meyer 1 12] i. We uke as granted that, with 
minor modifications, the results proved on the whole line can be transposed to 
the torus. In that section, homogeneous and inhomogeneous operators will be 
treated similarly. 


III.l General formulation 

The general formulation of the problem is 


(?) 


( And u € T\: ip such that 

Lu s / {S*\ 

nth f 6 fi- jp and l a contUat cctflieifnl elliptic operator of 
ardcr 5 . 


Standard variational approximation (P A. Raviart and J M Thomas [14]) 
loads us to look for the solution of a weak problem in so called Galerkin approx- 
imation spaces i; . where ? is a scale related to \' f with 

A natural choice for 1 1 is l, = V** ' where I’' belongs to a mult (resolution 
analysis of 7{o.:p of the type described above Indeed, we then have the following 
inequality guaranteed if the involved nmhireoolutton analysis b r-regular. 

V s < r. 3 r > 0 V/ 6 H’ l|/ - lit */i| 2 < C 2* ;f ||/!! w . (3$) 

where D r » . j < 0. stands for tbe orthogonal protection on V,” 

Then a standard Galerkin approximation writes 

Find Uy ^ Vf C T- such 

(p,) { n l>r /.n v t« r = n v ,/. (id) 


irhcrf n v 7 , j < 0 stands for the extension operator from 17 to 7\; ( : -v 

This approach lead* us tc replace L by the approximation FlrW-llv^. The 
corresponding numerical algorithm* are reduced tc linear system solvers ooce a 
basi« of V* ha< been chosen (P.A. Raviart et al [14)1 

tymbol 9 foe flit perkxLiatiuti on [0 t] ; We recall ih*\ dbmt ' * 

2 3} * 1/3 Am Wj hnf (kit s ipin’4^ m I). 


IS 


Numerically, an optimal choice for the expansiou basis is the wavelet basis 
of l p * : {♦? 0 .#£ . A € A£“ 1 }V since the corresponding stiffness matrix, is sparse 
and can be easily (i.e usiue diagonal matrices) uniformallv preconditioned (S 
Jaffard ($)). 


III. 2 A different approximation 


The main purpose of this paper is to define a different approximation of L and 


A;" - (defined in (31)) and defining if — 

spnn{£- *?<,. 

A € A,,}, we get 

an approximation of a as 


*,= £ + 

(11) 

as a; 



Indeed, 



« P = z - 1 n v */« 

% a 

(42) 


where Pp, is the projection or. Cf orthogonal to If. 

This formulation of u f shows that the convergence of u, towards u wlien 
p — oo is straightforward since the set t '* is a family of Galerkin spaces for the 
suitable spare of definition of m 

Moreover, the stability of the algorithm is a direct consequence from the 
classical preconditioning properties of wavelet base expansions (S. Jaffard [!?]) 


III.S General scheme 

Tlie numerical algorithm derived from the previous section is now presented 
in its collocation version We call /. the set of points t ZPl'O. ^l) 3 Then. 
J; = k € />} is the two dunensional regular grid of scale 2“ J Hated to 

10. 1 ? 

For the numerical implementation, we assume that the space V* is such that 
«ay coaftaaoas/aacftoa / £ I j* is intmki fuottsJf dtfmtd ly its rohet oa Me «rf 
rf points J t This assumption (satisfied by the even order spline multiresolution 
we Will use in the numerical tests) allows us to define the collocation projection. 
Cyp. from the net of real sequences t/i to VI” as 

C» ’((/*)♦!/,) = / <=>/€!,’ and f\2~ r k) = /». Vl e 1 T 

As soon as we define f 6 l'** by its coordinates on the basis the opera- 
tor C\t appears as a discrete convolution operator involving a so-called inter- 
polant filter /*,(£) Tlie operator C^i is also a convolution operator called the 

^A|«io. Aj, m o^„A ; »Hh. in the pniodiaed framenorfc. A, ■ ( \ * J"-’ (i - t ) .« £ 

» Ol°.> , f‘ , .» € r*;.. 


pom' value operator and involves the point value filter defined from 

Obviously. /»,(() and P\’i,(Z ) are inverse. Let us remark however, 
that the point value operator always exists as soon as the functions ♦** are 
continuous. 


For the implementation, we therefore replace D v * bv CY * in (42) and define 
therefore u T as 

«. * r l Cv*f 

Given the point values of f on the grid points J r . the algorithm provides the 
values of u P on the same grid. More precisely, the algorithm can be presented 
as follows: 

1. The input of the procedure is the set of values from which the 

interpolate functiou f y € 1JT is constructed using l+,i*Y 

h ~ 5Z f e*^r* 

u»» 

2. f r is then decomposed into the wavelet subspaces II’ ~ . 0 < j < p - 1 and 
to P » 

/ ~ •+• fe: 

»€V»-' 

when- coo — H\ »(/) 

P 

3* u r then become* 

*r = X + C U 

»e.vf-' 

aliere = (2*’I“ vx) Here. = r;o/<r((»| for non-homogeneous 
operators Tor homogeneous operator c^ should be given Note that in 
that case 0 ( 0 ) = 0 and / should have at least * vanishing moments: the 
fact that Cf 0 should be given corresponds to the ill posed property of the 
initial problem in L 2 

4. u r is then expanded in term? of the set of functions f?V k € l ? ) using 
the tree algorithms related to 0* and r* t as 

“r = 51 c r*V 

5. Finally, the grid poiut value? of a. on J, are estimated using the point 
value filter PV : if) 
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It appears dearly that various precalculations should be performed In the 
first step, the interpolant filter related to 4 r * must be known: for the second 
step, the orthofooal muhiresolution analysis quadrature mirror filters M, should 
be used and. for the fourthstep the corresponding btortbogona) multiresolution 
quadrature mirror filters P, are required; finally, the point value filter related 
to ry is used for the last step. 

To be more precise, we have to make some remarks that help to reduce the 
complexity and storage. For steps one and two. teiisorial properties can be used 
in a wry classical way to reduce the 2D- algorithms involved in 2(dim(l x 

lD-algorithms It is then enough to know the one-dimensional interpolant filter 
related to and the one dimensional quadrature mirror fillers m*. ( € E 

For steps four and five, where the filters P,. and the point value filter P\ ~ are 
involved, we can note apply this simplification and we have a full 2D*problem 
\Ye are now able to sununariie all there precalculations in the following step 

0 * 


0. The computation of the following filters is performed (this is presented for 
the spline muHuesoiution analysis case) 

-Interpolant filter related to >*,♦. I^ f : analytical formulae in on* di- 
mension are available in V Perrier and C Basdevant [IS]. 

-Orthogonal ID multircsolution analysis filters m , : analytical formu- 
las are also available 

•Filters P, of the BiorthogonaJ Multi resolution Analysis: There filters 
are constructed from to, and formula (23) and (24). In fact only P, 0 
and Pi, i have to be computed since we have P,i c.(f. {;) = Pc j.({ : ,« 4 ) 
Point value filter. P\'~ related to ryj This filter is computed from 

formula (25) and the analytical expression of .*?, c and ♦{£ ) We have 
successively 

l/'2») £ f r (u k : ** * . X 6 Bt 3 . 

£ ( £ hi* ’*■*'') -r.. € j, . (43: 

*€/, W M* / 

Practically 43 is truncated according to a prescribed precision. Thu is 
possible because decreases fast 

Remarks: 

• One should again emphaairt that the entire algorithm is based rn con 
volution operator* 1 hanks to the periodic boundary conditions, the 

convolutions can cither be performed directly or using a discrete Fourier 
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transform. The implementation presented in this paper uses the Fourier 
transforms since it is optimal for non compactly supported filters on non 
adapted spaces of approximation. 

III.4 Detailed Algorithm 

This section is devoted to the structure of the code. Basic tools, inch as Fast 
Fourier Transforms. Convolution/Decimaiion algorithms, or Term by terra mul- 
tiplications are not described. 

As can be seen from the general scheme presented above, the main code 
involves only two mote elaborate routines that will be called the /‘rece/ru/e; 
metiae (step 0). and the tree algorithm routines. The tree algorithm routines 
may or may not use the teoscrial structure. They will be called consequently !D 
Free Algorithm- D (steps 2) and W A'ea Ttntonnl Trxt Alfvntkm-I 
(step 4) where - D and -I stand for direct aad inverse. We recall that the steps 
1 and $ are convolutions and the renonnaluatioo performed in step 3 is term 
by term multiplication 

The tree algorithm routines are becoming very classical aad therefore we 
will not describe them either. Note however, that since only convolutions are 
performed in our algorithm, we only use the discrete Fourier transform of the 
wavelet coefficients (and not the corresponding values) at every scale, that re- 
duces significantly the complexity 

W> now give the detailed description of the main program (ELL IP ) and of 
the precalculus program (PASCAL) in pseudo code. 

The following example sketches the structure our programs 

[OUTPUTS] -Progran ( 2BPUTS ) 

* Cobb eats 

* Body oi Progran: 

< 

[OUTPUTSl]* Subprogram (IIPOTS) 

IIPUTS2 = OUTPUTS! 

[OUTPUTS] * Subprograms ( IIPUTS2) 

} 

Our variable descriptors beats some resemblance to the C language as well 
as to the MAT LAB conventions 


Til. 4.1 Preliminary computations 

The symbols •. and .1 used to present this program are borrowed from 
MAT LAB and mean respectively, the transposition, the matrix product the 
term by term product, and the term by term divUiou We also use the following 
u sub sampling operator an* defined as: If a is a 2D array of sire V > V 
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b^l : 2* : 2?. 1 : 2® : 2P) is ft new array of sire 2*~ n x y~ n given by ali./l r- 

aiTi.2.\ T j). (i.jl 6 l f . 

Prof rim PRECAL 

[QMFBIV.TAUTV] * PRECAL(p,paax,QNFV,PIIV,SV,CV, SIGMA) 

SIIPUT: 

*p •> index of the approximation space Vp is which the 

• elliptic problem is solve. 

•pat* -> index of the approximation space Vpaax in eluch the 

• precomputation of TAUTILDE is done (it depends on the 

• prescribed precision) . 

•QWV -> structure containing the quadrature airror filters in 
t one diaension: 

• QMFV.nO -> ID array containing the qaadratare airror filter 

• coefficients associated to the scaling Inactions; 

• sise(QMFV.aO)->2‘p; QMFV.aO(i) « u0(i/2*p), 

• i belong to <0, . . . ,2*p-l> . 

• QMFV.nl -> ID array containing the quadrature airror filter 

• coefficients associeted to the wavelet ; 

• eise(QMFY.nl)->2*p; 

• OSFV.nl(i) * al(i/2‘p) , i belong to {0 2‘p-l>. 

SPIIV -> ID array; sise<PlIV)->2*paax; where paax is given 

S and paax>p; PIIV(i) « the value of the Fourier transfora 

• of the ID scaling function at the point i, x belong to 

• {0, . . . ,2'paax-l}. Used to ceapute the value of tautilde on 

s the finer grid. 

*SV -> 2D array containing the sampling of the function 

• S used for biorthogoaal filters; 

• Sixe(SV)->(2*paax X 2‘paax); SV(i.j)«S(i/2*paax.j''2*paax). 

• (i.j) belong to <0 2 'paax- 1) "2. 

*CV -> 2D array containing tha sampling of the function 

• S(2v)/(2*s S(e)) Suse(CV) -> (2*p X 2*p); CV(i,J) » 

• S(2i/2*p,2j/2'p>/ S(i/2*p,j/2'p) . 

• <i.J) belong to <0. . . . ,2*p-l}*2. 

•SIGMA -> 2D array containing tha sampling of the tyabol 
» of tho operator Sase(SXCMA) -> (2*paar X 2‘paax); 

• SIGMA(i.j) ■ signa( i/2 *paax, j/2*p). 

• (i.j) belong to <0, . . . .2*paax-l)'2. 

•OUTPUT: 

SQKFBIV -> structure containiag the biorthogoaal quadratura 

• airror flltars related to the tautilde functions: 

• phFilW. PTILDE0-> 2D tab containing tho biorthogonal 

•JO 



* 

9 

* 

9 Computation of the filter* Ptilde 

• 

QIFBIV . PTILDEO » ((QHFV.nO) • * (QHFV.nO)) .♦ C V; 

QEFBIV. PTILDEX » ((QHFV.nO)* * (QKFV.srt)) ./ 

(2't SV(l:2'(p-pnax) :2'f*u, l:2‘(p-*pnax) :2'paax); 
QKFBIV PTILDE3 • ((QHFV.nl)* • (QHFV.nl)) ./ 

(2** SW(l:2*(p-pnn*):2*pnas.l:2*(p-pnax):2‘pnax); 

* 

* Computation «1 the point value filter related to TAUTILDE 

• 

Turm * ( (2*(p«) • P1IV* • PIIV) .• S$)./ SICIU; 

T1DTW * Per iodise (TAUTU , p); 

Thf subroutine Periodint is not described here, but it is a straight forward 
tranecnptKw of 43 

m.4.2 Mein Program 


quadrature eirror filters essocieted 
to the seeling functions; 
sise(QHFBIV.PTILDE0)->(2‘p X 2*p): 
QHFBIV. PTlLDEi-> 2D teb eenteiaisg the biortkegonel 
filters essociet *d to the first nr- 
•let ; sise( QHFBIV. PTILDEl )-> (2'p X 2‘p>. 
QHFBIV. PTILDE2-> sane ns QHFV. PTILDEl for the sseond 
wavelet (set conputod QKFBIV. PTXLDE2 * 
QHFV. PTILDEl transposed). 

QHFBIV. PTILDBl-> sene es QHFV. PTILDB1 for the third 

■evelet . 


.If# i« Prvfrvm Elhp 

CVX3 • Kllip(FX,QHFV,FIV.QBFBIV,T4irTV) 

•IIPOT: 

•FX -> 2D errey containing the sanplisg of the tenet ion 

• f; Sise(Fl) -> (2“p X 2*p); FX(i,j) ■ f (i/2*p. j/2*p) , 

9 (i,j) belong to (0, — ,2*p-l}“2. 

•FIV -> ID tab of data containing tho interpolation filter 
« related to PIX_pO, sise<FXV)->2*p; 

•QHFV -> see PBECAL 
•QHFBIV -> see PBECAL 
•TmV -> see PBECAL 

• 
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* 

•ODTWT: 

•UX ”> 2D ur»y cwtuiiaf tkt stapling of the approx i- 

• -nation u_p; tise(U])->(2'p X 2 p), OI(i,j) « 

• u(i/2'p,j/2*p) , (i.j) belong to <0 2*p-l)*2. 

PTitm? *> 2D array containing the values of the Fotirier 

• transior* of the scaling function TAUTHDE at level p 

•TQIP01ARY DATA : 

*FV -> 2D tab containing; tbo fft of FX; Sise(FY) -> 

« (2‘p I 2"p); 

«CPV -> 2D tab containing tbo fft of acaling coeffi-ient of 

• FX; Siae(CPV) -> <2*p X 2*p); 

•DJtf -> Structure of 2D array containing tho Four i ox transfon 
f of tho gavelet coefficients; tize(DJV) -> (2'p X 2*p); 
iCTILDEPW 

• -> sane as CPV for UX: 

•UV -> 2D tab containing tho fft of DX; 

i 


stop 0 

[FW] = Fast Fourier Transform(FX) 
step 1 

[CPV] * FV.eFIV (Ten by ton multiplication} 
stop 2 

[DJW] « 2D Tentorial Troe Algorithm.® (CPV.QHFH) 
step 3 

CDJV) * DJW. e(2* j*) (Term by term multiplication) 

stop 4 

tCTILDEPV) • 2D Ion Tensorinl Treo Algorithm.! (DJV.QHFBl) 
step S 
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[UV] « CTILSEPV . *TAUTM (Term by term mltiplication) 

CUX] * Inverse Fut Fourier Transfom(UV) 

) 

III.S Storage and Complexity Analysis 

As the computation is dearly separated into precalculations and actual imple- 
mentation of the algorithm, we will also present the storage and complexity 
aaalyses separating the two parts. One should remember that the preealeth 
lation is done once and for all while, as it will be the case in section IV, the 
algorithm can be applied iteratively. 

We will not discuss the complexity related to one-dimensional computations 
as well as the storage connected to one-dimensional arrays since both can be ne- 
glected in our bidimensional implementation All the evaluations are performed 
for .V = dim i; = y x 2* 

Storage 

Ptrmamt»( itonjt (pnco/rsiafiea* J The structures QNFBXV and TAUTtf rep- 
resent four bidimensional arrays of siae .V 

Ttmpornry siorvgt («ctunl The storage related to bidimensional 

arrays can be reduced to one arrays of siie X. 

F inally. the total memerp a serf corresponds to five arrays of sire «V 

Complexity analysis 

Prtcdlcmlts: The computation of the four arrays in the structure QXFBIV is 
done in C x .V. 

The computation of Pl~» is performed in C * X operations The value of 
C depends on the precision of the calculation 
,lfam program. 

Faat Four Aar Tran* for* and Inverse Fast Fourier Transform in- 
volve C x X log(.V) operations 

The complexity of the Tam by tarn multiplication is .V. 

Tiraa algoritlurt) and Trea algoritha-l are based on convolution and 
decimation operators. These procedures involve therefore C x .V opera- 
twos 

Therefore the fofaf complexity is 0( V k>gt ,V ) ) 

In the following section *r use these programs iteratively to solve the 2D 
Burgers equation after reducing it to a cascade of elliptic problems. We would 
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like to emphasize moreover, that our approach can he also used to solve equa- 
tions involving homogeneous pseudo-differential operators. A characteristic ex- 
ample is y-A u ~ f with periodi c bound ary conditions on [0, l) J . We have 
L — v'-A and therefore <r(() = y'ifjf + The most natural choice for 5 is 

5(f) — 2ysin } [fi/2! -f sin ! (^}/2) and ooe easily checks that the hypotheses of 
theorem II. 2 are then satisfied. The algorithms presented previously can be 
used (see Pj Ponenti [16]). 


IV Numerical Application: Resolution of the 
2D Burgers equations 


In this section, we will use the periodiaed Battle-ternaries multiresolution anal- 
ysis of splines of order m (see P.G. Lemarie [9]) The existence of collocation 
projectors related to the spline breakpoints requires splines of even order and 
the value m * ? will hr used in the applications 

As described in J. Ltandrat ct al [II], any parabolic equation of the type 


— i©» •+ f»(u) = 0 
t» = tie for t = (• 


144) 


0 < t < Tmar, 0 < z < 1 


where I? is a differential operator of even order with positive symbol and 

G is generally a nonlinear function of u and its derivatives, can be numerically 
approached using a classical finite difference time discretisation scheme followed 
by a variational approximation of the resulting elliptic problems. We show now 
that the approach developed in the previous sections can be used efficiently to 
provide this approximation 

Following J. Liaadrat et al 'll] w first introduce a segmentation (f, 
of [0. 7 mor] (i.e a sequence {t-.lJi- such that C = <© < »i <...<*«= Fmar) 
and now look for a sequence of functions of the t variable such 

that u <n, (r) is an approximation of «(*.{„). 

With At, = <»+ 1 — <»• 0 < a < Af . and considering first <441 as an ordinary 
differential equation in time, a standard finite-difference discretisation leads to 
the following iterative equation: 


c , U n+,) s /■( «'•’•• Al„ Af«_,.(*(u ,n ) (45) 

where C-. is a step forward operator that together with F is determined by 
the choice of the finite-difference approximation of the time-dependent ordinary 
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different i&l equation We always assume that this approximation is at least semi 
implicit for the linear part l-c and explicit for the nonlinear part. Therefore 
-o > 0 such that C n = (1 4 oAt n Le) where I stands for the identity operator. 
By hypothesis. eo(w) > O.Wu, - and then C„ has always a symbol hounded from 
below by l. 

Hence, assuming that {u 1 "” 4 *./ = 0 »} and {G(*' n ~ : ,1 = 0. ...» are 

known, the resolution of (45) falls under the scope of paragraph III. 5 and the 
resolution of (44) can be therefore performed iteratively. 

The bidimensiooal Burgers equation writes, with u = (u., t»j): 

-* Tu.w = t»Au 
v(0.II = utl.t) 

u = no for t = 0 

, 0 < I < Tmai 0 < i < 1. 

Choosing a constant step segmentation of O.Trnar] (ie., a segmentation 
such that 3At such that VO < n < M,t r . = an implicit Crank- Nicholson 
time scheme (second order) for the linear term (vAs). and an explicit second 
order Adams-Bashforth scheme for the nonlinear term (Vn u) we get 



and the solution can l>e written as 




I 


with 


511+1 




To fall completely under the scope of paragraph 11U. one should be able 
to evaluate the point values of the nonlinear term of (47) We used the sim- 
plest method available that comists, as classically done in spectral methods 
(C Canuto et al. [3.1. to “apply the nonlinear operator on the grid points’. 
More precisely the approximation of we used is PV(Vu u .v 0) ) = 

CyUCrf x (V u 'j Cv-etu 1 )) where x is a letm by term multiplication of finite 
sequences and Cy* is the collocation projection introduced in section 111.2. 

Then, for each time step n At. the problem clearly belongs tc> the class of 
elliptic problems studied in section (HI) with £ = / - vk A and / = 2u " ' - 

u 1 1 *). The iterative form of the equation induces 

some modifications of the general scheme presented in III 3 and we therefore 
provide the full scheme for an iteration of the Burgers approximation scheme. 


I 
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0. The inputs are the values of )l. ( 

1. F' lJ,) = J.' iJ,)- *,Pl(V„- VX-V) » twnputed u dncnM 

above. 

2. f^'irh the function of 1 r ~ interpolating the values F^yJ r ) is constructed 

* s/ T = E 

3. /*p is then decomposed iuto the wavelet subspaces H’ * . 0 < j < p - 1 and 
\ 7 as 

*o? -1 

where - n - = U v r(F?) 

4. ti/** 1 th*n beo?in« 

“r" + " = £ 2--'.:r;.vr 5:-ra: 


where t»J = (y/.-f,)’ 


5. is then expanded in terms of tk 1 set cf functions {r* t . k £ ! r ) v 


M r ~ 2- f ?' it 

t«x. 


We also get = E* e /, f i*^V* 

and s E» t i/r*£?,V 

6. From the point values of r*\J F ). fir- \J r \ and we compute 

(u' n ' , (;,)) and its gradient using the corresponding point value filters. 

7. Finally using the values of <«.,< J P )i and of it* gradient, we get the values 
of (u 1 ’ 1 * 1 ( J r )| = (u ’ ,+ *(</,) - «- R • J F ' i and its gradient 

IV. 1 Storage and Complexity Analysis 

As described above the numerical code implement* the elLptic solver in ar. it- 
erative process. Since the time step St is constant, the characteristics of the 
elliptic solver does not depend on the time index n. Then, the soNer precalcu- 
lations related to the whole parabolic problem are the same as the ones related 
to the elementary elliptic solver t.aee section 111 a). This applies to the storage 
and to the complexity as well. 



As it has l >mi shown in the previous section however, extra work, not con- 
nected to the elliptic solver itself but to the computation of the right hand side 
term of the iterative equation (47) is required. This extra work h related to the 
storage of the fields at the different time steps involved in the three level time 
step Adams-Ba&hforth Crauk- Nicholson scheme and to the point value evalu- 
ation of the derivatives involved in the nonlinear part. Again, it can be split 
into permanent and temporary storage as well as in precalculation and main 
program extra work 

Storage (in addition to the elliptic solver storage) 

Terms* sat sloragt (pttcolemUUons): Ore extra structure containing the 
point values must be stored in one bidimensronal array of site A*; the 

structure containing the point values (/,) t$ given by transposition of the 
previous one. 

Temporary storm ft (actual at fonthm; Thetwofields tr n >. ti*"- 1 and u 
can He handled using three arrays of site A\ 

Complexity analysis (in addition to the elliptic solver complexity) 

PrtcolcmlmUom: The computation of is performed in C x A’ opera- 

tions where, as in section HI 5 the value of C depends on the precision of the 
calculation. 

.Vim program: The addition of complexity 1 $ related to steps (I). (5). (6) 
and (7). Since these steps involve convolutions and term by term products, the 
added complexity is again CX log V 

Finally, the total winery ased is 7 arrays of sire X 
The total compftntp is QUC logtAJ) 

Obviously^ the total complexity of the whole resolution is M times the com- 
plexity of one time step resolution 

TV.2 Numerical Results 

Test case on an translation invariant problem 

The validation of our code ha* been performed on an translation invari- 
ant problem constructed using for the initial condition (ur,. uc,)(r : . x 3 ) s 
(sin(2ir ]), 0). Indeed with such an initial condition, the solution remains r* 
translation invariant. 

For an easy comparison to the well documented paper of C Basdevant et al 
:f we used u = 10~ 2 /r 

As explained in C. Bnsdevant ei al (1], the pertinent quantities are 
. 8u fa 

ms^sujH\ r=sup* 0 (0.5. f): 
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and tmj> defined such that 


}^-(0.5,f,„)| = ms 

Table 4S exhibits the numerical results obtain using various values for thr 
time step At The maximum time step numerically acceptable was At = O.OC75. 
In each case, the values of ms are computed by interpolation and the correspond* 
ing values of t-«, are deduced. The comparison with the expected theoretical 
values (first column) show* that oar method competes favorably with the ma- 
jority of the schemes presented in C. Basdevaut et al [1] A complete study of 
the time step sire dependence of the results connected to the stability analysis 
of the parabolic algorithm will be published later. 


Kiad 

A? 

0.0005 

oooi 

0.(025 

0.005 

-304.0103 

ms 

-304.630? 

-305.727 

-309.4354 

—316.5434 

0. $55237 


0.233 

0.252 

0.23 

0 243 


Test case on a first diagonal translation invariant problem 

Our second test case is performed on a first diagonal translation invariant prob- 
lem constructed using (u^,.ae J )tr.,X}) = (sin(2r(ri j- xj).sin(2r(r ; •+ xjt) 
Again, the solution can be compared to the reference solution of C- Basdevnrit 
et ai. [1] thanks to a 45* rotation and to a time dilation of factor 2. However, 
according to onr reference axes, it is obviously a fully bidimensional solution. 

Figures 1.2 and 3 show the isoline values of the numerical approximation* 
computed with 6t = .001 at f = 0. < = 0.15 and t = 0.50. The fitst diagonal 
translation invariance is kept and we obtain the values ms = 249 052$ and t^, = 
0 123. The expected theoretical values are -304.0103 for ms and 0.1276185 for 
If.?. This is not as good as before but one should note that the resc4ut>on in 
the direction perpendicular to the front axis is now half the one in our previous 
calculations. 


Since the ultimate application of all this work is the development of adaptive 
algorithms (i.e. the development of algorithms handling approximation spaces 
of reduced dimension adapted to the solution regularity (see for instance Pj. 
Ponenti (ICji.we have estimated, at various times, the wavelet basis adapted 
to the approximation and defined as the lowest cardinal m = $ spline wavelet 
basis preserving the £ 3 norm of the approximation with a precision of 10~*. 
The <c4unins of table (49) show foe each Kale 0 < J < 7 the number of wavelet* 
selected in tlie adapted basis related to the approximated solution at various 
times. It appears that, compared to the full basis of IV <la#t column), these 
bases have a drastically reduced cardinal (we defined the rate of compression 
u — ~) e w n if the gradients of the solution fib up • large 

domain made of two complete lines of the plane (sec figures 1. 2 and 3) 






Figure 3: Approximated solution, t, = 0 50 
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Conclusion 

In this paper, we have proposed an inversion scheme for elliptic problems 
based on biorihogonal wavelets The approximation of elliptic problem solutions 
is constructed and leads to stable and fast numerical algorithms. 

Numerical tests related to the approximation of the parabolic Burgers equa- 
tions transformed into a cascade of elliptic problems are provided. 

The approximation scheme is baaed on convolution operators and can there- 
fore be theoretically used in the framework of adapted space* of approximation 
As mentioned however, the nice Censorial product structure that enforces nu- 
merical efficiency is fragile and is generally destroyed when applying the scheme 
directly Other approximations for the step forward operator, should allow 
one to use efficiently this approximation in a general context of adapted multi- 
dimension spaces of approximation. 
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